http://adni.loni.usc.edu/data-samples/access-data/
Image taken from http://www.martinos.org/neurorecovery/technology.htm
cort_thickness in extrantsr packagelabel roi 1 X1002 “left caudal anterior cingulate” 2 X1003 “left caudal middle frontal” 3 X1005 “left cuneus” 4 X1008 “left inferior parietal” 5 X1017 “left paracentral” 6 X1022 “left postcentral” 7 X1023 “left posterior cingulate” 8 X1024 “left precentral” 9 X1025 “left precuneus” 10 X1027 “left rostral middle frontal” 11 X1028 “left superior frontal” 12 X1029 “left superior parietal” 13 X1030 “left superior temporal” 14 X1031 “left supramarginal” 15 X1035 “left insula”
source('utils.R')
source('combat.R')
modelData = read.csv('modelData.csv')
head(modelData)
subject age sex dx scanner 1 002_S_0413 76.4329 F Normal 3_PHILIPS_Intera_2 2 002_S_0559 79.4658 M Normal 3_PHILIPS_Intera_2 3 002_S_0729 65.2986 F MCI 3_PHILIPS_Intera_2 4 002_S_0816 70.9534 M AD 3_PHILIPS_Intera_2 5 002_S_0954 69.5041 F MCI 3_PHILIPS_Intera_2 6 002_S_1018 70.8055 F AD 3_PHILIPS_Intera_2
mod = model.matrix(~age+factor(sex)+factor(dx), data=modelData)
ctData = read.csv('imageData.csv')
head(ctData)[,1:10]
subject X1002 X1003 X1005 X1008 X1017 X1022
1 002_S_0413 2.349515 1.957340 1.6771022 2.929418 1.893073 1.733991
2 002_S_0559 2.814481 1.768518 0.9173002 2.297948 1.612640 2.071636
3 002_S_0729 2.788202 2.108902 1.6343228 3.154604 2.219242 1.984391
4 002_S_0816 2.753477 2.168249 1.7955870 2.880065 1.973897 2.042263
5 002_S_0954 2.523273 1.664635 1.0189855 2.077749 1.236037 1.468252
6 002_S_1018 2.596688 2.216399 1.9206076 2.424231 1.744081 1.800574
X1023 X1024 X1025
1 2.365209 1.968747 2.474627
2 2.606214 1.838712 2.345573
3 2.578613 1.782150 2.055732
4 2.188748 1.990839 2.992091
5 1.329828 1.468782 1.623318
6 2.165901 1.947617 2.611285
img = t(ctData[,-1]) head(img)[,1:10]
[,1] [,2] [,3] [,4] [,5] [,6] [,7]
X1002 2.349515 2.8144805 2.788202 2.753477 2.523273 2.596688 1.833424
X1003 1.957340 1.7685180 2.108902 2.168249 1.664635 2.216399 1.696198
X1005 1.677102 0.9173002 1.634323 1.795587 1.018986 1.920608 1.228951
X1008 2.929418 2.2979484 3.154604 2.880065 2.077749 2.424231 2.712218
X1017 1.893073 1.6126404 2.219242 1.973897 1.236037 1.744081 1.760396
X1022 1.733991 2.0716363 1.984391 2.042263 1.468252 1.800574 1.712189
[,8] [,9] [,10]
X1002 2.960386 3.191757 1.6549346
X1003 2.220978 2.294089 1.1845660
X1005 1.638375 1.760091 0.6901224
X1008 2.727885 2.280574 1.2541459
X1017 2.105533 1.890207 1.0815922
X1022 2.095173 1.480126 0.8755437
harmonized = combat(dat=img, batch=modelData$scanner, mod=mod)
[combat] Performing ComBat with empirical Bayes [combat] Found 14 batches [combat] Adjusting for 4 covariate(s) or covariate level(s) [combat] Standardizing Data across features [combat] Fitting L/S model and finding priors [combat] Finding parametric adjustments [combat] Adjusting the Data
combat() function returns a list.head(harmonized$dat.combat)[,1:10]
[,1] [,2] [,3] [,4] [,5] [,6] [,7]
X1002 2.075969 2.5403388 2.514211 2.479313 2.2494910 2.322785 1.5601570
X1003 1.604147 1.3928650 1.774367 1.829854 1.2947712 1.890831 1.3247258
X1005 1.251037 0.6043039 1.215196 1.363586 0.6807333 1.459647 0.8726338
X1008 2.242658 1.6391016 2.457074 2.191413 1.4269869 1.755756 2.0304396
X1017 1.414838 1.1241720 1.749303 1.497199 0.7381837 1.263199 1.2756135
X1022 1.172603 1.5470296 1.459725 1.521212 0.8741849 1.260857 1.1429860
[,8] [,9] [,10]
X1002 2.686305 3.134451 1.7216344
X1003 1.886106 2.388280 1.3679891
X1005 1.219925 1.922309 0.8532241
X1008 2.053663 2.445597 1.4071359
X1017 1.632919 2.067534 1.2650911
X1022 1.579223 1.739706 1.1057018
harmonized$gamma.hat versus harmonized$gamma.starharmonized$delta.hat versus harmonized$delta.starmodelData$sex = factor(modelData$sex) modelData$dx = factor(modelData$dx) modelData$scanner = factor(modelData$scanner)
preComBat = left_join(modelData, ctData, by='subject') preRInsula = lm(X2035 ~ age + sex + dx + scanner, data=preComBat) summary(aov(preRInsula))
Df Sum Sq Mean Sq F value Pr(>F) age 1 4.222 4.222 24.157 2.48e-06 *** sex 1 0.297 0.297 1.701 0.19427 dx 2 2.007 1.004 5.743 0.00402 ** scanner 13 12.456 0.958 5.482 6.07e-08 *** Residuals 138 24.119 0.175 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
harmonizedData = as.data.frame(t(harmonized$dat.combat)) harmonizedData$subject = ctData$subject
postComBat = left_join(modelData, harmonizedData, by='subject') postRInsula = lm(X2035 ~ age + sex + dx + scanner, data=postComBat) summary(aov(postRInsula))
Df Sum Sq Mean Sq F value Pr(>F) age 1 1.625 1.6250 10.377 0.00159 ** sex 1 0.550 0.5497 3.511 0.06309 . dx 2 1.799 0.8994 5.744 0.00402 ** scanner 13 0.529 0.0407 0.260 0.99566 Residuals 138 21.610 0.1566 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Das, Sandhitsu R, Brian B Avants, Murray Grossman, and James C Gee. 2009. “Registration Based Cortical Thickness Measurement.” Neuroimage 45 (3). Elsevier: 867–79.
Fortin, Jean-Philippe, Nicholas Cullen, Yvette I Sheline, Warren D Taylor, Irem Aselcioglu, Philip A Cook, Phil Adams, et al. 2018. “Harmonization of Cortical Thickness Measurements Across Scanners and Sites.” NeuroImage 167. Elsevier: 104–20.
Fortin, Jean-Philippe, Drew Parker, Birkan Tunc, Takanori Watanabe, Mark A Elliott, Kosha Ruparel, David R Roalf, et al. 2017. “Harmonization of Multi-Site Diffusion Tensor Imaging Data.” Neuroimage 161. Elsevier: 149–70.
Johnson, W Evan, Cheng Li, and Ariel Rabinovic. 2007. “Adjusting Batch Effects in Microarray Expression Data Using Empirical Bayes Methods.” Biostatistics 8 (1). Oxford University Press: 118–27.